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Abstract: Rapeseed is considered as one of the most important oilseed crops in the world. Vegetable 
oil obtained from rapeseed is a valuable raw material for the food and energy industry as well 
as for industrial applications. Compared to other vegetable oils, it has a lower concentration of 
saturated fatty acids (5%-10%), a higher content of monounsaturated fatty acids (44%-75%), and 
a moderate content of alpha-linolenic acid (9%-13%). Overall, rapeseed is grown in all continents 
on an industrial scale, so there is a growing need to predict yield before harvest. A combination 
of quantitative and qualitative data were used in this work in order to build three independent 
prediction models, on the basis of which yield simulations were carried out. Empirical data collected 
during field tests carried out in 2008-2015 were used to build three models, QQWR15_4, QQWR31_5, 
and QQWR30_6. Each model was composed of a different number of independent variables, ranging 
from 21 to 27. The lowest MAPE (mean absolute percentage error) yield prediction error corresponded 
to QQWR31_5, it was 6.88%, and the coefficient of determination R 2 was 0.69. As a result of the 
sensitivity analysis of the neural network, the most important independent variable influencing the 
final rapeseed yield was indicated, and for all the analyzed models it was "The kind of sowing date 
in the previous year" (KSD_PY). 


Keywords: winter rapeseed; artificial neural networks; machine learning; yield prediction; yield 
simulation; MLP (multi-layer perceptron) network; sensitivity analysis; precision agriculture; 
sustainable agriculture; sustainability 


1. Introduction 

Vegetable oil obtained from rapeseed is an important element of food and a balanced diet. It is 

also used for the production of biodiesel, in the form of rapeseed methyl esters (RME) and is also 
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an essential industrial raw material. These main applications of rapeseed are associated with a high 
content of protein, carbohydrates, and triacylglycerol (TAG) [1-4]. In view of the above, the cultivated 
land used for rapeseed must compete between food, industrial, and energy purposes [5,6]. Thus, 
production planning and forecasting of rapeseed yield is crucial in the whole cultivation process. This is 
directly related to the financial policy and sales deadline strategy adopted by rapeseed producers. 

The 2017 world rapeseed production amounted to 76,238,340 tons, with an average yield of 
2.19 t ha -1 on the total crop land of 34,740,403 ha. The largest rapeseed producers in the world include 
Canada (21,328,000 tons), China (13,274,000 tons), India (7,917,000 tons), France (5,200,000 tons), 
Australia (4,313,230 tons), Germany (4,275,600 tons), and Poland (2,697,265 tons) [7]. 

There are many factors that influence the size and quality of rapeseed and other crops. These include 
climatic and soil factors, sowing-cultivation-harvesting technologies, plant variety, and fertilization 
level [8]. Moreover, the mentioned factors are directly or indirectly related to each other, and their 
effect on the final crop yield is often not thoroughly studied. Therefore, prediction models supporting 
decision-making processes throughout the entire plant cultivation stages are increasingly used for 
data analysis to determine the potential yield before harvest. Such actions are aimed at optimizing 
production in terms of consumption of all production means which in turn maximize yield and financial 
profits while minimizing inputs. 

In recent years, an increase in the use of prediction models in agriculture has been observed. 
In particular, tools based on artificial intelligence—artificial neural networks (ANN), giving significantly 
lower predichon errors than in case of statistical methods—are very popular. Therefore, crop yield 
models are implemented in computer applications for precision agriculture and are becoming an 
important element of decision support systems [9-18]. 

Most often, quantitative data is used for model building, which makes it easier to choose the 
modeling method. However, there are situations that require a different forecasting approach, which 
complicates the model building process. The use of combined data (quantitative and qualitative) may 
in consequence give a better result than in case of only one type of data (only quantitative or only 
qualitative). Therefore, the authors of this work, encouraged by previous research results [11,19,20], 
have undertaken the construction and analysis of winter rapeseed yield models based on quantitative 
and qualitative data. The produced models can be used to simulate the crop in the current agrotechnical 
season before harvesting. The main assumption is to achieve MAPE (mean absolute percentage error) 
of up to 7%. 

Some of the developed forecasting models predict yields too late, only a few days before 
harvesting [21]. Al-Gaadi et al. [21] attempted to correlate the final potato tuber yield, i.e., basing on 
samples taken three days before harvest, with the results of vegetative indices during plant growth. 
The results of the study showed that the difference between the predicted yield values and actual, 
observed values was small, as it did not exceed 15%. Of course, crop yield control from air and satellite 
level is valuable in agricultural production [22], but requires the use of specialized equipment and 
software. In addition, due to the frequent association of yield with local spectral characteristics of 
plants, such models can only be used in the region for which they were developed [23]. 

Yield predictions based on readily available agronomic and meteorological data made at relatively 
early and key stages of plant growth and development have a greater contribution to conscious 
plantation control [20]. During the winter rapeseed vegetation, five critical stages of seed yield 
formation are indicated, i.e., autumn vegetation, spring re-vegetation, formation of the main and side 
shoots, flowering, growth of pods, and maturation. While yield prediction in the first two phases may 
seem too early, the estimation of seed yield after formation of the main and side shoots, flowering, 
growth of pods and maturation stage, provides interesting, preliminary information about the actual 
potential yielding of plants in situ. 

In the construction of forecasting models, not all available independent variables determining the 
scale of plant yielding should be taken into account. A big practical problem is combining the vegetation 
process of plants according to the widely used BBCH scale (Biologische Bundesanstalt, Bundessortenamt 
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und CHemische Industrie) with the dynamically changing production conditions [24-26]. Operating 
on the BBCH scale allows for relatively reliable analyses, however precise collection of such data 
during vegetation is often simply impossible. This would require farmers to regularly observe the 
plantations and meticulously record the progress of vegetation. Only those farms that manage a small 
area and cultivate one or two varieties would have a chance to succeed. In the case of cultivating 
several varieties on a dozen or several dozen areas of fields simultaneously, very accurate monitoring of 
development phases causes problems. Therefore, in many works, the authors of the models use specific 
dates of the calendar year, which correspond to key developmental phases in typical years [10,11,20]. 

2. Materials and Methods 

The research place and material was described in [19]. However, additional qualitative data 
were used to build predictive models (Table 1, Table 2). All field data took into consideration 8 years 
of research (from 2008 to 2015), from winter rapeseed fields located in Poland, in the central and 
southwestern part of Greater Poland (Figure 1). In total, data from 328 fields were collected, which 
were divided into two sets (Table 1). Set A (292 fields) from the first seven years was the basis for 
model development, while set B (36 fields) from the last year of research was used to validate the 
constructed models and determine yield prediction errors. 




Figure 1. Research area — Greater Poland, Poland (redraw from [19]). 

Table 1. The number of productive fields of winter rapeseed divided into two sets, A and B [19]. 






Set A 




Set B 

Year 

2008 

2009 

2010 

2011 

2012 

2013 

2014 

2015 

Number of fields 

32 

49 

48 

50 

45 

28 

40 

36 


Basic meteorological data was necessary for the construction of predictive models, and the data 
was obtained from meteorological stations located near fields from which empirical data was collected. 
Air temperature and atmospheric precipitation were divided into those obtained in the previous year 
(from September 1 to December 31) and from the current agrotechnical year. Meteorological data for 
the current year have been additionally divided into shorter periods depending on the model: from 
January 1 to April 15, from April 1 to April 30, from May 1 to May 31 and from June 1 to June 30. 

The neural predictive model construction has been made on the basis of three predicted dates 
of a calendar year: April 15, May 31 and June 30. They have been named respectively QQWR15_4 
(QuantitativeQualitativeWinterRapeseedl5_April), QQWR31_5 and QQWR30_6. The prediction dates 
proposed allow easier control of winter rapeseed yield. They were chosen approximately to the course of 
key vegetation phases in average years: i.e., April 15—inflorescence emergence. May 31—development 
and formation of pods, June 30—technical maturity. 

The models took into consideration those factors (independent variables) that affect crop yields and 
are widely available to agricultural producers (Table 2). The following approach for winter rapeseed 
yield prediction allows forecasting and simulation of estimated yields directly before harvesting, in the 
same agricultural year. 
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Table 2. Data structure in the neural prediction models. 


Symbol 


Variable Name 




The Scope of Data 

Measure 

QQWR15_4 

QQWR31_5 

QQWR30_6 

Quantitative data 

R9-12 PY 

mm 

The sum of precipitation from September 




63-234 


1 to December 31 of the previous year 






The average air temperature from 





T9-12_PY 

°C 

September 1 to December 31 of the 

+ 

+ 

+ 

4.9-9.4 



previous year 





Rl-4 CY 

mm 

The sum of precipitation from January 1 




59-185 


to April 15 of the current year 




Tl-4 CY 

°c 

The average air temperature from 




-0.4-4.9 



January 1 to April 15 of the current year 




R4 CY 

mm 

The sum of precipitation from April 1 to 




8.7-60.4 


April 30 of the current year 




T4 CY 

°c 

The average air temperature from April 




5.9-12.2 



1 to April 30 of the current year 




R5 CY 

mm 

The sum of precipitation from May 1 to 




14.2-132.5 


May 31 of the current year 




T5 CY 

°C 

The average air temperature from May 1 




11.8-16.2 



to May 31 of the current year 




R 6 CY 

mm 

Total precipitation from June 1 to June 30 




15-121 


of the current year 




T 6 CY 

°c 

The average air temperature from June 1 




14.2-19.6 



to June 30 of the current year 




N_LY 

kg ha -1 





0-41 

the previous year 

+ 

+ 

+ 

N CY 

kg ha -1 

The sum of N fertilization - spring in the 

+ 

+ 

+ 

0-175 


current year 





P205 CY 

kg ha -1 

The sum of P 2 O 5 fertilization in the 

+ 

+ 

+ 

0-104 


current year 





K2CLCY 

kg ha -1 

The sum of K 2 O fertilization in the 




0-234 

current year 




MGO CY 

kg ha -1 

The sum of MgO fertilization in the 

+ 

+ 

+ 

0-298 


current year 





S03 CY 

kg ha -1 





14-115 


+ 

+ 

+ 


current year 





B CY 

gha - 1 

The sum of B fertilization in the current 

+ 

+ 

+ 

0-3.66 


year 





CU_CY 

MN CY 

gha - 1 

The sum of Cu fertilization in the current 
year 

The sum of Mn fertilization in the 

+ 

+ 

+ 

10-487 

70-600 

gha - 1 

+ 

+ 

+ 


current year 





MO_CY 

gha - 1 

The sum of Mo fertilization in the 




0-60 

current year 




ZN CY 

gha - 1 

The sum of Zn fertilization in the current 

+ 

+ 

+ 

9-226 


year 





DOF PY 

tha - 1 

The dose of organic fertilizer in the 

+ 

+ 

+ 

0-41 



previous year 





Qualitative data 

KOF PY 

word 

The kind of organic fertilizer in the 

+ 

+ 

+ 

No, 

Livestock manure. 



previous year 




Chicken manure 

LF_PY 

word 

The liming of a field in the previous year 

+ 

+ 

+ 

Yes, 

No 







Aggregate cultivator. 

ST_PY 

word 

The soil tillage in the previous year 

+ 

+ 

+ 

Shallow tillage. 

Deep tillage. 







Ploughing 

KSD PY 

word 

The kind of sowing date in the 




Early, 


previous year 




Optimal Tardy 

Winter rapeseed. 








Winter wheat. 

TF_PY 

word 

The type of forecrop in the previous year 

+ 

+ 

+ 

Winter triticale. 

Winter barley. 

Spring barley. 







Oat, Rye 


"+"—the variable exists in the model, —the variable does not exist in the model. 


Tables 1 and 2 present independent variables in qualitative and quantitative research that were the 
basis for the construction of each of the three neural models. The choice of the network topology was 
guided by the results of earlier studies [10,11,20]. In addition, the network's ability to approximate and 
generalize, considering the network error rate, was taken into account. Each constructed neural model 
QQWR15_4, QQWR31_5, and QQWR30_6 was a result of learning 10,000 networks, one of which was 
selected with the best parameters. Basic measures of network quality were determined—standard 
deviation, mean error, deviation error, mean absolute error, quotient of deviation, and correlation. 
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The best model was chosen based on the lowest mean absolute error value and the largest correlation 
value. Automatic network designer from Statistica v7.1 (StatSoft Inc., Tulsa, OK, USA) was used for 
neural network building. The 292 data set was divided into three sets: training, validation and test. 
The sets were divided into 70%-15%-15% proportion (204-44-44 cases in each set). 

In the next stage, the obtained neural models QQWR15_4, QQWR31_5, and QQWR30_6 were 
validated, comparing the results of model prediction with real data from set B (36 cases from 2015). 
Data from set B were not considered in construction of QQWR15_4, QQWR31_5, and QQWR30_6 
neural models. In order to determine the prediction errors of winter rapeseed yields, the difference 
between the observed and the predicted value was performed. This method is widely described in 
the literature [20,27-30]. Determination of the quality of predictions was based on the error measures 
(Equations (l)-(4)): 

RAE—relative approximation error; 


RAE = 


rUy-hY 

\ LUiyt) 2 


RMS—root mean square error; 


RMS - 


L'Uy-yi) 2 


MAE—mean absolute error; 


MAE = li,\y<-. 


i=l 


MAPE—mean absolute percentage error; 


MAPE = - 
n 


i n 

ty 

n 

i=l 


Vi ~ Vi 


Vi 


x 100% 


( 1 ) 


( 2 ) 

(3) 

(4) 


where 

n —number of observations, 

y,—actual values obtained during research, 

y,—values given by the model. 

The last stage of the research was the sensitivity analysis of the neural network. This analysis gives 
the opportunity to examine all independent variables in each produced neural model QQWR15_4, 
QQWR31_5, and QQWR30_6. As a result of this test, the information about the importance (rank) 
of each measure is obtained. The greater the numerical value of a given feature (error quotient), 
the greater the impact on winter rapeseed yield. All values below 1 indicate a low impact on the yield 
of the given independent variable and can be removed from the model. 

3. Results 

According to the presented methodology of neural model building, 1 of 10,000 models was 
selected for analysis considering a specific forecast date. These models were chosen based on the best 
qualitative parameters of the neural network. Detailed information on the quality and structure of 
selected neural models QQWR15_4, QQWR31_5, and QQWR30_6 are presented in Table 3. 










Agronomy 2019 , 9, 781 


6 of 11 


Table 3. The quality and structure of the neural models produced. 



QQWR15_4 

QQWR31_5 

QQWR30_6 

Neural network 

MLP 

MLP 

MLP 

structure 

21:33-15-9-1:1 

25:37-16-16-1:1 

27:39-10-1:1 

Learning error 

0.0515 

0.0914 

0.0932 

Validation error 

0.0898 

0.0987 

0.0837 

Test error 

0.1033 

0.1042 

0.1386 

Mean 

3.3627 

3.3626 

3.3626 

Standard deviation 

1.0703 

1.0703 

1.0703 

Average error 

-0.0126 

0.0924 

0.0205 

Deviation error 

0.4087 

0.5578 

0.5548 

Mean Absolute error 

0.2870 

0.4353 

0.4118 

Quotient deviations 

0.3818 

0.5211 

0.5183 

Correlation 

0.9245 

0.8534 

0.8552 


MLP—multi-layer perceptron. 


The next stage of the analysis was validation of the produced models. The observed yield from 
set B was compared to the results of prediction of all produced models. For this purpose, formulas 
(1-4) were applied and the results of model validation are included in Table 4. 


Table 4. Measures prediction of analyzed neural models. 


Model 

RAE (-) 

RMS (t ha 1 ) 

MAE (t ha 1 ) 

MAPE (%) 

QQWR15_4 

0.098 

0.43 

0.37 

9.87 

QQWR31_5 

0.068 

0.29 

0.25 

6.88 

QQWR30_6 

0.076 

0.32 

0.28 

7.69 


Netzuork Sensitivity Analysis 

The sensitivity analysis of the produced neural models was the last stage of calculations. For each 
model QQWR15_4, QQWR31_5, and QQWR30_6 the error quotient and the rank were determined, 
they are presented in Table 5. 


Table 5. Sensitivity analysis of neural networks. 


Variable 



Model 





QQWR15.4 

QQWR31_5 

QQWR30_6 


quotient 

rank 

quotient 

rank 

quotient 

rank 

R9-12_PY 

1.549 

4 

1.077 

5 

1.029 

10 

T9-12_PY 

1.224 

9 

1.007 

20 

1.027 

11 

R1-4_CY 

1.314 

5 

1.028 

9 

1.018 

14 

T1-4_CY 

1.129 

13 

1.027 

10 

1.012 

20 

R4_CY 

- 

- 

1.041 

7 

1.036 

8 

T4_CY 

- 

- 

1.012 

17 

1.021 

13 

R5_CY 

- 

- 

1.068 

6 

1.078 

6 

T5_CY 

- 

- 

1.237 

2 

1.228 

2 

R6_CY 

- 

- 

- 

- 

1.016 

17 

T6_CY 

- 

- 

- 

- 

1.012 

19 

N_LY 

1.099 

15 

1.016 

15 

1.010 

21 

N_CY 

1.107 

14 

1.007 

19 

1.000 

26 

P205_CY 

1.027 

18 

1.013 

16 

0.989 

27 

K2CLCY 

1.012 

21 

0.994 

24 

1.003 

25 

MGOCY 

1.232 

7 

1.009 

18 

1.017 

16 

S03_CY 

1.081 

16 

1.038 

8 

1.018 

15 

B_CY 

1.150 

11 

1.027 

11 

1.008 

22 

CU_CY 

1.014 

20 

1.005 

21 

1.004 

23 

MN_CY 

1.018 

19 

0.994 

25 

1.015 

18 

MO_CY 

1.637 

3 

1.016 

14 

1.130 

5 

ZN_CY 

1.078 

17 

0.999 

23 

1.003 

24 

DOF_PY 

1.144 

12 

1.002 

22 

1.029 

9 

KOF_PY 

1.228 

8 

1.026 

12 

1.040 

7 

LF_PY 

1.169 

10 

1.016 

13 

1.025 

12 

ST_PY 

1.294 

6 

1.098 

4 

1.144 

4 

KSD_PY 

2.273 

1 

1.251 

1 

1.435 

1 

TF_PY 

1.846 

2 

1.204 

3 

1.150 

3 
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By considering the results obtained in the previous stages, a decision was made to perform 
additional analyses and visualizations between the observed yield values and the yield values 
predicted by the QQWR15_4, QQWR31_5, and QQWR30_6 models. The results of the analyses are 
presented in Figure 2. The QQWR31_5 model obtained the highest coefficient of determination, 
the value of which was R 2 = 0.69. Slightly lower R 2 values were calculated for models QQWR15_4 and 
QQWR30_6, which were 0.47 and 0.60, respectively. 


Yield observed vs Yield predicted by model QQWR15_4 
y = 1.25 + 0.71*x; p = 0.00; R 2 = 0.47 



(a) 

Yield observed vs Yield predicted by model QQWR31_5 
y = 0.98 + 0.75*x; p = 0.00; R 2 = 0.69 



(b) 

Figure 2. Cont. 
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Yield observed vs Yield predicted by model QQWR30_6 
y = 1.32 + 0.66‘x; p = 0.00; R 2 = 0.60 



(c) 

Figure 2. Relation between the observed and predicted yield (t ha -1 ) by model: (a) QQWR15_4, 

(b) QQWR31_5 and (c) QQWR30_6. 

4. Discussion 

Many different modeling methods are applied to predict crop yield. Regardless of the method 
chosen, it is necessary to collect empirical data from the past. It is assumed that the models will 
have a lower forecast error as the number of years and the number of data on which they are based 
increase [10]. A large number of factors affecting crop yields determine the choice of the optimal method 
of analysis. Many of the independent variables can only be obtained from specialized tests. Apart 
from data availability, another problem in modeling is their format and structure. While quantitative 
data are most often used, available qualitative data are often difficult to implement in produced 
models. Therefore, the decision was made to integrate quantitative and qualitative data. Each model 
produced as part of the following research work, i.e., QQWR15_4, QQWR31_5, and QQWR30_6 is 
based on both types of data and allows making yield forecasts in the current agrotechnical year before 
harvesting. In addition, it was assumed that a good model will have a low MAPE and will be up to 
7%. Artificial neural networks, which are non-linear modeling methods, were used to implement the 
assumed prediction error. As a rule, they give better forecast results than classical mathematical models 
and enable combining quantitative and qualitative data. Most often, MLP (multi-layer perceptron) 
topology network is used in forecasting, which was implemented to build three models QQWR15_4, 
QQWR31_5, and QQWR30. Each model is based on 21, 25, and 27 independent variables, respectively, 
which include 5 qualitative variables (Table 2). 

Validation of produced neural models for predicting winter rapeseed yields is based on determining 
the quality of forecasts using four measures, i.e., RAE, RMS, MAE, and MAPE. The error values for 
all produced neural models are shown in Table 4. However, the most commonly used indicator for 
determining the forecast error is MAPE. The lowest error values of this indicator were obtained only 
for the QQWR31_5 neural model based on the MLP network with the structure 25: 37-16-16-1: 1, which 
was 6.88%. The same time period (January-May) obtained the lowest error in prediction of winter 
wheat yield in work [20]. In the following work, the MAPE error rate was assumed to be up to 7%. 
This means that the other two models QQWR15_4 and QQWR30 did not achieve the assumed error 
rate because their MAPE was higher and amounted to 9.87 and 7.69, respectively. 

The network sensitivity analysis was another element of performed calculations. It took into 
account all neural networks produced. The independent variable "The kind of sowing date in the 
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previous year" (KSD_PY) (Table 5) obtained the highest rank 1 in each network, which is consistent 
with literature sources [31-34]. In the QQWR31_5 model, the error quotient was 1.2510. Interestingly 
enough, the independent variable, "The average air temperature from May 1 to May 31 of the current 
year" (T5_CY) came second in this model, which took a lower value of 1.2367. A similar situation 
was also observed in QQWR30 model. As it was investigated in [20], the dominant influence of 
atmospheric factors in predicting winter wheat yield did not correspond completely to prediction of 
winter rapeseed yield. In the QQWR15_4 model, the independent feature "The type of forecrop in 
the previous year" (TF_PY) came second. However, as indicated earlier, QQWR15_4 and QQWR30 
did not achieve the assumed MAPE rate. Despite the differences in the results of individual models, 
it can be explicitly indicated that qualitative data had a dominant impact on the forecast of winter 
rapeseed yield. Other authors studying this species [35-37] have also confirmed that sowing date is a 
very important feature. Rapeseed is very sensitive to date of sowing, because its dynamic development 
takes place at a temperature exceeding 12 °C. In the area studied, such temperatures usually occur at 
the turn of the second and third decade of August. Sowing that is too late, taking place in September 
or October, causes lower aboveground matter development in plants. Their growth cone contains 
insufficient buds of future generative organs. On the other hand, too early sowing poses the risk of 
excessive outgrowth before winter, which is characterized by a large above-ground matter, susceptible 
to damage and freezing. Failure to meet the proper sowing date is almost always associated with 
poorer wintering of plants. The maximum delay in sowing, which does not influence reduction in 
yield, is 5-7 days. 

Comparing the predicted yield to observed yield (Figure 2), the best fit was for QQWR31_5 model, 
with R 2 =0.69. This confirms the right assumption of the maximum MAPE of up to 7%. 

5. Conclusions 

The results of the present work clearly indicate the possibility of using neural models to predict 
the yield of winter rapeseed in practice, basing on quantitative and qualitative data. For this purpose, 
particularly useful is the QQWR31_5 model, whose MAPE is 6.88% and its determination coefficient 
R 2 is 0.69. It is the only model that met the assumptions of the work regarding the maximum MAPE 
of up to 7%. In addition, neural network sensitivity analysis revealed that "The kind of sowing date 
in the previous year" (KSD_PY) had the greatest impact on winter rapeseed yield. Such result was 
obtained for all analyzed models. The neural models presented in the paper can be applied in precision 
agriculture as a key factor in decision making systems. 

Further research should be aimed at optimizing the produced neural models in two directions. 
As part of the first direction, independent variables should be reduced while simultaneously maintaining 
the MAPE rate. The second direction is to obtain more empirical data from regions with different 
soil and climatic conditions. Such actions will aid in expanding the field scope of model practical 
application and should positively affect the accuracy of forecast. 
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